%% 
% 
% Turn detection code for 
% "Vibrio cholerae motility in aquatic and mucus-mimicking environments"
% by
% Marianne Grognot, Anisha Mittal, Mattia A. Mah’moud, Katja M. Taute*
% Rowland Institute at Harvard
% 100 Edwin H Land Blvd
% Cambridge, MA 02142
% USA
% *taute@rowland.harvard.edu
% 
% The article is accepted at Applied and Environmental Microbiology.
% An earlier version of the manuscript is available on bioRxiv at
% https://doi.org/10.1101/2021.07.06.451398


% BUGTURNS 
%
% INPUT :
%   V : Speed structure with (filtered) trajectories to analyze
% OUTPUT :
%   T : Turn analysis structure, for each bug in V. Empty is not entering
%   analysis
%       T(bug#,1) is a logical matrice or runs (each run has a colum, with
%       true values at corresponding lines/index in the speed structure of V)
%       T(bug#,2) gathers the index of run start(s)
%       T(bug#,3) gathers the index of run end(s)
%
function T=BugTurns(V)

% Taking info from input structure
Speeds=V.Speeds;
fps=V.Parameters.fps;
medV=cell2mat(cellfun(@(x) (nanmedian(x(:,9))), Speeds, 'UniformOutput', 0));

% Threshold to enter analysis
SpeedThreshold=20;% average bacterial speed minimal to be included (um/s)
TimeThreshold=1;% number of second minimal to be included (s)
T=cell(size(Speeds,1), 3);

disp(strcat('ANALYSIS BUGTUMBLE. Minimal Time :',num2str(TimeThreshold),' s. Minimal Speed :',num2str(SpeedThreshold),' um/s'))

for i=1:length(Speeds)

    if ( size(Speeds{i},1)>(TimeThreshold*fps) ) && ( medV(i)>SpeedThreshold )
        disp(['Analyzing Bug ' num2str(i)]);
        
        X=Speeds{i}(:,2:4);% trajectory from Speeds is the filtered one, already scaled
        dX=Speeds{i}(:,6:8);
        dAlpha=Speeds{i}(:,end);        

        [R,runstart,runend]=RunTurn_Vcholerae(X,dX,dAlpha,fps);
        
        %       
        T{i,1}=R;
        T{i,2}=runstart;
        T{i,3}=runend;
    end
end


function [R, runstart, runend]=RunTurn_Vcholerae(X, dX, dAlpha,fps)
%
ThreshFactorA=8;%
medA=nanmedian(dAlpha);

loop=1;
while loop<4
    loop=loop+1;  

ThreshAngle=ThreshFactorA*medA;
n=size(X,1);

%initialize run and turn matrices. Each column is an event. Overestimate
%number of events, then later throw out the rest.
R=false(n,round(n));

%declare beginning to be run.
R(1:2,1)=logical( [1;1]);

%keep track of tripped  thresholds
T=logical(zeros(2,3));

%logical index saying whether dAlpha is larger than a threshold angle
%for angles
Lalpha=dAlpha>ThreshAngle;

%collect start times of runs and tumbles
runstart=[];
runend=[];


jr=1; %run counter
jt=0;%tumble counter
i=3;
while i<n-3
    %if currently running, check for turn (i.e. the end of a run)
    if R(i-1,jr)==1

        %tumble (run end) criterion
        if Lalpha(i)
            %declare new turn: turn counter is increased by one
            jt=jt+1;
            runend=[ runend i];
        else
            %otherwise declare continued run
            R(i,jr)=1;
        end
        i=i+1;
    else
        %if you're turning, check whether the beginning of a run
        %should be called
        
        if ~any(Lalpha(i:i+0))
            %start new run
            jr=jr+1;
            R(i,jr)=1;%
            runstart=[runstart i];
            i=i+1;
        else
            i=i+1;
            
        end
        
        
    end
    
end
R(:,sum(R)==0)=[];

% Now that a first iteration defined the "Runs", compute medA on them only
dAlpha_Runs=[];
for j=1:size(R,2)
    dAlpha_Runs=[dAlpha_Runs; dAlpha(R(:,j))];
end

medA=nanmedian(dAlpha_Runs);% New medA for next iteration
end



